import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import make_classification, load_iris

# 设置随机种子以确保结果可重现
np.random.seed(42)

##hhfai
def pca_manual(X, n_components=2):
    """
    手动实现PCA算法
    """
    # 1. 标准化数据（中心化）s
    X_standardized = StandardScaler().fit_transform(X)

    # 2. 计算协方差矩阵
    cov_matrix = np.cov(X_standardized.T)

    # 3. 计算特征值和特征向量
    eigenvalues, eigenvectors = np.linalg.eig(cov_matrix)

    # 4. 对特征值排序并选择前n_components个特征向量
    sorted_indices = np.argsort(eigenvalues)[::-1]
    eigenvectors_sorted = eigenvectors[:, sorted_indices]
    eigenvectors_subset = eigenvectors_sorted[:, :n_components]

    # 5. 将数据投影到新的特征空间
    X_pca = X_standardized.dot(eigenvectors_subset)

    return X_pca, eigenvalues, eigenvectors


def pca_sklearn(X, n_components=2):
    """
    使用sklearn实现PCA
    """
    # 标准化数据
    X_standardized = StandardScaler().fit_transform(X)

    # 创建PCA对象
    pca = PCA(n_components=n_components)

    # 拟合数据并进行转换
    X_pca = pca.fit_transform(X_standardized)

    return X_pca, pca


# 生成示例数据
def generate_sample_data():
    """生成一个高维分类数据集"""
    X, y = make_classification(
        n_samples=300,
        n_features=10,  # 原始特征维度
        n_informative=8,
        n_redundant=2,
        n_clusters_per_class=1,
        random_state=42
    )
    return X, y


# 可视化结果
def plot_results(X_pca, y, title="PCA降维结果"):
    """可视化降维后的数据"""
    plt.figure(figsize=(10, 8))
    scatter = plt.scatter(X_pca[:, 0], X_pca[:, 1], c=y, cmap='viridis', alpha=0.7)
    plt.colorbar(scatter)
    plt.xlabel('主成分 1')
    plt.ylabel('主成分 2')
    plt.title(title)
    plt.grid(True, alpha=0.3)
    plt.show()


def main():
    # 方法1：使用生成的数据
    print("=== 使用生成的高维数据 ===")
    X, y = generate_sample_data()
    print(f"原始数据形状: {X.shape}")

    # 使用手动实现的PCA
    X_pca_manual, eigenvalues, eigenvectors = pca_manual(X, n_components=2)
    print(f"手动PCA降维后形状: {X_pca_manual.shape}")
    print(f"特征值: {eigenvalues[:5]}")  # 显示前5个特征值

    # 使用sklearn的PCA
    X_pca_sklearn, pca = pca_sklearn(X, n_components=2)
    print(f"Sklearn PCA降维后形状: {X_pca_sklearn.shape}")
    print(f"解释方差比: {pca.explained_variance_ratio_}")
    print(f"累计解释方差: {np.cumsum(pca.explained_variance_ratio_)}")

    # 可视化结果
    plot_results(X_pca_manual, y, "手动PCA实现")
    plot_results(X_pca_sklearn, y, "Sklearn PCA实现")

    # 方法2：使用真实数据集（鸢尾花数据集）
    print("\n=== 使用鸢尾花数据集 ===")
    iris = load_iris()
    X_iris, y_iris = iris.data, iris.target
    print(f"鸢尾花数据形状: {X_iris.shape}")

    # 使用sklearn PCA
    X_iris_pca, pca_iris = pca_sklearn(X_iris, n_components=2)
    print(f"降维后形状: {X_iris_pca.shape}")
    print(f"解释方差比: {pca_iris.explained_variance_ratio_}")
    print(f"累计解释方差: {np.cumsum(pca_iris.explained_variance_ratio_)}")

    # 可视化鸢尾花数据
    plot_results(X_iris_pca, y_iris, "鸢尾花数据集PCA降维")

    # 绘制解释方差比（碎石图）
    def plot_explained_variance(X, title="解释方差比"):
        """绘制解释方差比的碎石图"""
        pca_full = PCA().fit(StandardScaler().fit_transform(X))
        explained_variance = pca_full.explained_variance_ratio_

        plt.figure(figsize=(10, 6))
        plt.plot(range(1, len(explained_variance) + 1),
                 np.cumsum(explained_variance), 'bo-', alpha=0.8)
        plt.xlabel('主成分数量')
        plt.ylabel('累计解释方差比')
        plt.title(title)
        plt.grid(True, alpha=0.3)
        plt.axhline(y=0.95, color='r', linestyle='--', alpha=0.7, label='95%方差')
        plt.legend()
        plt.show()

    plot_explained_variance(X, "生成数据的解释方差比")
    plot_explained_variance(X_iris, "鸢尾花数据的解释方差比")


if __name__ == "__main__":
    main()